(ns curve-fitting
  (:use (incanter core stats charts))
  (:use (util random)))

;;; This is for the curve fitting example in Bishop's book
;;; I've implemented linear polynomial regression with regularization
;;; Results match the book

(defn parse-data []
  (let [lines (.split (slurp "data/curve-fitting.dat") "\n")
        fields (map #(.split % " ") lines)]
    {:xs (map #(Double/parseDouble %) (map first fields))
     :ts (map #(Double/parseDouble %) (map second fields))}))

(def *data* (parse-data))
(def *xs* (:xs *data*))
(def *ts* (:ts *data*))

(defn ** [x n] (reduce * (repeat n x)))
(defn times2 [x] (+ x x))
(defn pow-list [xs i] (map #(** % i) xs))
(defn sum-pow-list [xs i] (reduce + (pow-list xs i)))

(defn regression-polynomial-lhs
  ([xs M]
     "generate a matrix with the following form
   | sum(pow(x,0)) sum(pow(x,1))  .. sum(pow(x,M))  |
   | sum(pow(x,1)) sum(pow(x,2))  .. sum(pow(x,M+1) |
   | ...           ...            .. ...            |
   | sum(pow(x,M)) sum(pow(x,M+1) .. sum(pow(x,2M)  |"
     (let [ff (fn [i] (sum-pow-list xs i))
           xs (map ff (range (inc (times2 M))))]
       (loop [xs xs n 0 rs []]
         (if (= n (+ M 1))
           (matrix rs)
           (recur (rest xs) (inc n) (conj rs (take (inc M) xs)))))))
  ([xs M lambda]
     "linear polynomial regression with regularization"
     (plus (regression-polynomial-lhs xs M)
            (mult (identity-matrix (inc M)) lambda))))

(defn regression-polynomial-rhs
  ([xs ts M]
     "generate the column vector in the following form:
  | sum(t)              |
  | sum(t * x)          |
  | ...                 |
  | sum(t * pow(x, M+1) |"
     (let [x-pow (fn [i] (reduce + (mult ts (map #(** % i) xs))))]
       (matrix (map (fn [i] (x-pow i)) (range (inc M)))))))

(defn regression-polynomial-coefficients
  ([xs ts M lambda]
     (mmult (solve (regression-polynomial-lhs xs M lambda))
            (regression-polynomial-rhs xs ts M)))
  ([xs ts M]
     (regression-polynomial-coefficients xs ts M 0)))


(defn polynomial [ws x]
  (let [xs (map #(** x %) (range (count ws)))]
    (reduce + (mult ws xs))))

(defn polynomial-function [ws]
  (fn [x]
    (polynomial ws x)))

(defn polynomial-regression
  "returns a function"
  ([xs ts M lambda]
     (let [ws (regression-polynomial-coefficients xs ts M lambda)]
       (fn [x] (polynomial ws x))))
  ([xs ts M]
     (polynomial-regression xs ts M 0)))

(defn polynomial-regression-plot
  ([M xs ts lambda]
     (let [f (polynomial-regression *xs* *ts* M lambda)
           plot (function-plot f 0 1)]
       (doall
        (map (fn [i] (add-pointer
                     plot
                     (nth xs i)
                     (nth ts i))) (range (count xs))))
       plot))
  ([M] (polynomial-regression-plot M *xs* *ts* )))

(defn signal-generator
  ([f m s]
     (fn [x]
       (let [g (gaussian-generator m s)]
         (+ (f x) (g)))))
  ([m s]
     (signal-generator (fn [x] (Math/sin (* 2 Math/PI x))) m s))
  ([f]
     (signal-generator f 0 1))
  ([]
     (signal-generator (fn [x] (Math/sin (* 2 Math/PI x))) 0 0.3)))



